In this report I will go through the use of Modern Portfolio Theory to optimize a portfolio of any number of stock holdings by measuring the annual returns against the annualized risk, while recording the Sharpe Ratio of each portfolio.

Select tickers that make up your portfolio and assign the price data.

tick <- c('SCHG','RSP','SPLG','IHI','TSM',
          'MSFT',   'EFX','FCX','BBJP', 'XLRE', 'DIS',  'PYPL',
          'AMZN',   'PSA',  'GOOG', 'VPU')

price_data <- tq_get(tick, from = '2020-01-01',
                     get = 'stock.prices')
## change the from date if you want further historical data, I like to do 
## around 10 years

## enter the weights of your portfolio make sure the weights correspond to the 
## order of tickers entered in tick above 
my_wts <- c(0.1463, 0.1269, 0.0661, 0.0646, 0.0518, 0.0456, 
            0.0408, 0.0407, 0.0397, 0.0387, 0.0358,
            0.0345, 0.0332, 0.0328, 0.0307, 0.0261)

Next transform the data using log scale

log_ret_tidy <- price_data %>%
  group_by(symbol) %>%
  tq_transmute(select = adjusted,
               mutate_fun = periodReturn,
               period = 'daily',
               col_rename = 'ret',
               type = 'log')

head(log_ret_tidy)
## # A tibble: 6 × 3
## # Groups:   symbol [1]
##   symbol date            ret
##   <chr>  <date>        <dbl>
## 1 SCHG   2020-01-02  0      
## 2 SCHG   2020-01-03 -0.00554
## 3 SCHG   2020-01-06  0.00628
## 4 SCHG   2020-01-07 -0.00138
## 5 SCHG   2020-01-08  0.00699
## 6 SCHG   2020-01-09  0.00914
log_ret_xts <- log_ret_tidy %>%
  spread(symbol, value = ret) %>%
  tk_xts()
## Warning: Non-numeric columns being dropped: date
## Using column `date` for date_var.
head(log_ret_xts)
##                    AMZN          BBJP           DIS          EFX          FCX
## 2020-01-02  0.000000000  0.0000000000  0.0000000000  0.000000000  0.000000000
## 2020-01-03 -0.012213314 -0.0105434629 -0.0115372638  0.003365160 -0.030771493
## 2020-01-06  0.014775871  0.0036622412 -0.0058189870  0.004260731  0.003898638
## 2020-01-07  0.002089436  0.0016233857  0.0003432508 -0.001883466  0.015444163
## 2020-01-08 -0.007839287  0.0004054758 -0.0020611691  0.021210391  0.009912407
## 2020-01-09  0.004787695  0.0060617598 -0.0039278697  0.011554637 -0.016832816
##                     GOOG          IHI         MSFT           PSA         PYPL
## 2020-01-02  0.0000000000  0.000000000  0.000000000  0.0000000000  0.000000000
## 2020-01-03 -0.0049193493 -0.007216658 -0.012529993  0.0107264290 -0.018131770
## 2020-01-06  0.0243581503  0.008345515  0.002581339  0.0035500870  0.012880971
## 2020-01-07 -0.0006242443 -0.002032792 -0.009159304 -0.0039712338 -0.004548771
## 2020-01-08  0.0078495086  0.004211534  0.015802594  0.0018241000  0.019414593
## 2020-01-09  0.0109839096  0.008816560  0.012415846  0.0004203794  0.006684815
##                      RSP         SCHG         SPLG          TSM           VPU
## 2020-01-02  0.0000000000  0.000000000  0.000000000  0.000000000  0.0000000000
## 2020-01-03 -0.0057035680 -0.005538329 -0.006312541 -0.033533995  0.0002833025
## 2020-01-06  0.0012991377  0.006281609  0.003161252 -0.011606873  0.0019817689
## 2020-01-07 -0.0008657874 -0.001380641 -0.002633594  0.016074890 -0.0025484543
## 2020-01-08  0.0029408196  0.006990197  0.005260586  0.007346080  0.0000000000
## 2020-01-09  0.0045672302  0.009140058  0.007318301  0.008137101  0.0045260057
##                     XLRE
## 2020-01-02  0.0000000000
## 2020-01-03  0.0072987978
## 2020-01-06  0.0002599757
## 2020-01-07 -0.0109661283
## 2020-01-08  0.0049755020
## 2020-01-09 -0.0002612056
tail(log_ret_xts)
##                    AMZN         BBJP          DIS           EFX          FCX
## 2023-05-25 -0.015102768  0.005059203 -0.010496122 -0.0029406852  0.012567550
## 2023-05-26  0.043475866  0.003224506  0.001700409  0.0168038732  0.033336402
## 2023-05-30  0.012822304 -0.006256986 -0.005337600  0.0005695070 -0.015652520
## 2023-05-31 -0.008916851 -0.008948578  0.001592894 -0.0102056274  0.003208421
## 2023-06-01  0.017999211  0.023820268  0.007136788 -0.0003355566  0.032376443
## 2023-06-02  0.011983006  0.019947305  0.024309854  0.0241602982  0.047352568
##                    GOOG           IHI         MSFT           PSA          PYPL
## 2023-05-25  0.022034300 -0.0088155387  0.037736809 -0.0009886295 -0.0401143218
## 2023-05-26  0.008647678 -0.0005653933  0.021160155  0.0125660382  0.0142155130
## 2023-05-30 -0.006318259 -0.0069989680 -0.005059560  0.0076790852  0.0299359560
## 2023-05-31 -0.010241585  0.0043564714 -0.008550618 -0.0195742226 -0.0009673906
## 2023-06-01  0.008073023  0.0118367012  0.012678421  0.0004587895  0.0169549415
## 2023-06-02  0.006891058  0.0091103701  0.008443439  0.0173486283  0.0143298233
##                      RSP         SCHG          SPLG          TSM          VPU
## 2023-05-25 -0.0006429918  0.020643492  0.0084440798  0.113372272 -0.013992142
## 2023-05-26  0.0083271410  0.019941832  0.0132425596  0.022140424  0.000569165
## 2023-05-30 -0.0017734904  0.002695155 -0.0004048676 -0.011988990 -0.003776673
## 2023-05-31 -0.0087716429 -0.004401194 -0.0050746075 -0.033806949  0.008884524
## 2023-06-01  0.0074928315  0.013567182  0.0089142089  0.002532545 -0.006602615
## 2023-06-02  0.0215908529  0.009639000  0.0150137936  0.001011286  0.011050602
##                    XLRE
## 2023-05-25 0.0022695565
## 2023-05-26 0.0115510465
## 2023-05-30 0.0030765122
## 2023-05-31 0.0064022358
## 2023-06-01 0.0002773844
## 2023-06-02 0.0208627421

Next calculate the mean returns, and covariance matrix

mean_ret <- colMeans(log_ret_xts, na.rm = TRUE)
print(round(mean_ret, 5))
##     AMZN     BBJP      DIS      EFX      FCX     GOOG      IHI     MSFT 
##  0.00031  0.00010 -0.00057  0.00050  0.00124  0.00070  0.00024  0.00089 
##      PSA     PYPL      RSP     SCHG     SPLG      TSM      VPU     XLRE 
##  0.00052 -0.00064  0.00032  0.00051  0.00038  0.00066  0.00012  0.00008
cov_mat <- cov(log_ret_xts, use= 'complete') * 252

print(round(cov_mat,4))
##        AMZN   BBJP    DIS    EFX    FCX   GOOG    IHI   MSFT    PSA   PYPL
## AMZN 0.1505 0.0363 0.0663 0.0614 0.0770 0.0907 0.0511 0.0921 0.0344 0.1106
## BBJP 0.0363 0.0395 0.0431 0.0384 0.0678 0.0405 0.0333 0.0411 0.0234 0.0467
## DIS  0.0663 0.0431 0.1393 0.0650 0.1175 0.0694 0.0539 0.0679 0.0382 0.0915
## EFX  0.0614 0.0384 0.0650 0.1292 0.0797 0.0641 0.0565 0.0675 0.0488 0.0861
## FCX  0.0770 0.0678 0.1175 0.0797 0.3267 0.0901 0.0763 0.0884 0.0489 0.1153
## GOOG 0.0907 0.0405 0.0694 0.0641 0.0901 0.1185 0.0575 0.0945 0.0386 0.0966
## IHI  0.0511 0.0333 0.0539 0.0565 0.0763 0.0575 0.0627 0.0617 0.0388 0.0755
## MSFT 0.0921 0.0411 0.0679 0.0675 0.0884 0.0945 0.0617 0.1164 0.0458 0.1036
## PSA  0.0344 0.0234 0.0382 0.0488 0.0489 0.0386 0.0388 0.0458 0.0816 0.0448
## PYPL 0.1106 0.0467 0.0915 0.0861 0.1153 0.0966 0.0755 0.1036 0.0448 0.2440
## RSP  0.0498 0.0402 0.0702 0.0617 0.1027 0.0596 0.0530 0.0613 0.0425 0.0719
## SCHG 0.0852 0.0413 0.0698 0.0675 0.0952 0.0833 0.0605 0.0886 0.0433 0.1041
## SPLG 0.0613 0.0379 0.0646 0.0590 0.0913 0.0660 0.0528 0.0699 0.0408 0.0790
## TSM  0.0725 0.0446 0.0617 0.0565 0.1147 0.0744 0.0515 0.0772 0.0284 0.0939
## VPU  0.0310 0.0293 0.0453 0.0511 0.0582 0.0407 0.0444 0.0461 0.0488 0.0455
## XLRE 0.0483 0.0365 0.0589 0.0642 0.0805 0.0552 0.0542 0.0613 0.0626 0.0669
##         RSP   SCHG   SPLG    TSM    VPU   XLRE
## AMZN 0.0498 0.0852 0.0613 0.0725 0.0310 0.0483
## BBJP 0.0402 0.0413 0.0379 0.0446 0.0293 0.0365
## DIS  0.0702 0.0698 0.0646 0.0617 0.0453 0.0589
## EFX  0.0617 0.0675 0.0590 0.0565 0.0511 0.0642
## FCX  0.1027 0.0952 0.0913 0.1147 0.0582 0.0805
## GOOG 0.0596 0.0833 0.0660 0.0744 0.0407 0.0552
## IHI  0.0530 0.0605 0.0528 0.0515 0.0444 0.0542
## MSFT 0.0613 0.0886 0.0699 0.0772 0.0461 0.0613
## PSA  0.0425 0.0433 0.0408 0.0284 0.0488 0.0626
## PYPL 0.0719 0.1041 0.0790 0.0939 0.0455 0.0669
## RSP  0.0661 0.0625 0.0600 0.0591 0.0504 0.0620
## SCHG 0.0625 0.0831 0.0668 0.0757 0.0450 0.0606
## SPLG 0.0600 0.0668 0.0595 0.0615 0.0463 0.0574
## TSM  0.0591 0.0757 0.0615 0.1553 0.0313 0.0506
## VPU  0.0504 0.0450 0.0463 0.0313 0.0669 0.0597
## XLRE 0.0620 0.0606 0.0574 0.0506 0.0597 0.0794

Now we will find some random weights for the holdings in the portfolio, as well as finding that weighted portfolio’s returns, risk, and Sharpe ratio.

wts <- runif(n = length(tick))

wts <- wts/sum(wts)
print(wts)
##  [1] 0.05115565 0.01187139 0.04695456 0.09506853 0.13586896 0.01062504
##  [7] 0.10156061 0.10008125 0.03874295 0.04424691 0.01679664 0.02361168
## [13] 0.12790694 0.10435312 0.06155783 0.02959793
sum(wts)
## [1] 1
port_returns <- (sum(wts * mean_ret) + 1)^252 - 1
port_returns # we see an ___ return
## [1] 0.1246907
port_risk <- sqrt(t(wts) %*% (cov_mat %*% wts))
print(port_risk) # we see a ____ 
##          [,1]
## [1,] 0.273082
# Since Risk free rate is 0% 

sharpe_ratio <- port_returns/port_risk
print(sharpe_ratio) # sharpe ratio is ____
##           [,1]
## [1,] 0.4566052

Once we have made sure the work is correct, we begin the function to run the 5000 (can use different number) trials. Firstly, create the blank vectors within the DF to store numbers.

num_port <- 15000

# Creating a matrix to store the weights

all_wts <- matrix(nrow = num_port,
                  ncol = length(tick))

# Creating an empty vector to store
# Portfolio returns

port_returns <- vector('numeric', length = num_port)

# Creating an empty vector to store
# Portfolio Standard deviation

port_risk <- vector('numeric', length = num_port)

# Creating an empty vector to store
# Portfolio Sharpe Ratio

sharpe_ratio <- vector('numeric', length = num_port)

Now we begin the loop

for (i in seq_along(port_returns)) {
  
  wts <- runif(length(tick))
  wts <- wts/sum(wts)
  
  # Storing weight in the matrix
  all_wts[i,] <- wts
  
  # Portfolio returns
  
  port_ret <- sum(wts * mean_ret)
  port_ret <- ((port_ret + 1)^252) - 1
  
  # Storing Portfolio Returns values
  port_returns[i] <- port_ret
  
  
  # Creating and storing portfolio risk
  port_sd <- sqrt(t(wts) %*% (cov_mat  %*% wts))
  port_risk[i] <- port_sd
  
  # Creating and storing Portfolio Sharpe Ratios
  # Assuming 0% Risk free rate
  
  sr <- port_ret/port_sd
  sharpe_ratio[i] <- sr
  
}

Now we store the values in a table

# Storing the values in the table
portfolio_values <- tibble(Return = port_returns,
                           Risk = port_risk,
                           SharpeRatio = sharpe_ratio)


# Converting matrix to a tibble and changing column names
all_wts <- tk_tbl(all_wts)
## Warning in tk_tbl.data.frame(as.data.frame(data), preserve_index, rename_index,
## : Warning: No index to preserve. Object otherwise converted to tibble
## successfully.
colnames(all_wts) <- colnames(log_ret_xts)

# Combing all the values together
portfolio_values <- tk_tbl(cbind(all_wts, portfolio_values))
## Warning in tk_tbl.data.frame(cbind(all_wts, portfolio_values)): Warning: No
## index to preserve. Object otherwise converted to tibble successfully.
head(portfolio_values)
## # A tibble: 6 × 19
##     AMZN   BBJP     DIS    EFX    FCX   GOOG     IHI   MSFT     PSA    PYPL
##    <dbl>  <dbl>   <dbl>  <dbl>  <dbl>  <dbl>   <dbl>  <dbl>   <dbl>   <dbl>
## 1 0.0165 0.110  0.104   0.0358 0.103  0.0304 0.00935 0.0899 0.0791  0.0146 
## 2 0.0906 0.0401 0.0984  0.0675 0.0482 0.0426 0.0663  0.0846 0.0365  0.102  
## 3 0.0116 0.112  0.142   0.0351 0.0758 0.103  0.0245  0.0942 0.0974  0.0354 
## 4 0.0968 0.0159 0.0820  0.107  0.0735 0.0781 0.0789  0.0743 0.0414  0.00941
## 5 0.0459 0.0504 0.128   0.107  0.0290 0.131  0.0271  0.109  0.00607 0.0223 
## 6 0.0291 0.0549 0.00565 0.0979 0.0537 0.0659 0.0539  0.110  0.0229  0.115  
## # ℹ 9 more variables: RSP <dbl>, SCHG <dbl>, SPLG <dbl>, TSM <dbl>, VPU <dbl>,
## #   XLRE <dbl>, Return <dbl>, Risk <dbl>, SharpeRatio <dbl>

Now we will utilize our weights inputted at the top and view their statistics.

my_wts <- my_wts/sum(my_wts)
print(my_wts)
##  [1] 0.17125132 0.14854267 0.07737329 0.07561746 0.06063444 0.05337703
##  [7] 0.04775840 0.04764134 0.04647079 0.04530025 0.04190565 0.04038394
## [13] 0.03886223 0.03839401 0.03593585 0.03055133
sum(my_wts)
## [1] 1
my_port_returns <- (sum(my_wts * mean_ret) + 1)^252 - 1
my_port_returns # we see our return 
## [1] 0.08020782
my_port_risk <- sqrt(t(my_wts) %*% (cov_mat %*% my_wts))
print(my_port_risk) # we see our risk
##          [,1]
## [1,] 0.254466
# Since Risk free rate is 0% 

my_sharpe_ratio <- my_port_returns/my_port_risk
print(my_sharpe_ratio) # sharpe ratio 
##           [,1]
## [1,] 0.3152006
my_info <- c(my_wts, my_port_returns, my_port_risk, my_sharpe_ratio)
my_labels <- c(tick, "Return", "Risk", "SharpeRatio")
my_port <- data.frame(my_labels, my_info)

Then find the portfolios with min variance, optimal Sharpe ratio, and max returns

min_var <- portfolio_values[which.min(portfolio_values$Risk),]
max_sr <- portfolio_values[which.max(portfolio_values$SharpeRatio),]
max_ret <- portfolio_values[which.max(portfolio_values$Return),]

Then can create some graphs and charts to see our data

Histograms that reflect holdings in the three selected portfolios

# min variance
p <- min_var %>%
  gather(tick, key = Asset,
         value = Weights) %>%
  mutate(Asset = as.factor(Asset)) %>%
  ggplot(aes(x = fct_reorder(Asset,Weights), y = Weights, fill = Asset)) +
  geom_bar(stat = 'identity') +
  theme_minimal() +
  labs(x = 'Assets', y = 'Weights', title = "Minimum Variance Portfolio Weights") +
  scale_y_continuous(labels = scales::percent) 
## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
##   # Was:
##   data %>% select(tick)
## 
##   # Now:
##   data %>% select(all_of(tick))
## 
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
ggplotly(p)
### max returns
p <- max_ret %>%
  gather(tick, key = Asset,
         value = Weights) %>%
  mutate(Asset = as.factor(Asset)) %>%
  ggplot(aes(x = fct_reorder(Asset,Weights), y = Weights, fill = Asset)) +
  geom_bar(stat = 'identity') +
  theme_minimal() +
  labs(x = 'Assets', y = 'Weights', title = "Maximum Returns Portfolio Weights") +
  scale_y_continuous(labels = scales::percent) 

ggplotly(p)
### optmizied sharpe ratio
p <- max_sr %>%
  gather(tick, key = Asset,
         value = Weights) %>%
  mutate(Asset = as.factor(Asset)) %>%
  ggplot(aes(x = fct_reorder(Asset,Weights), y = Weights, fill = Asset)) +
  geom_bar(stat = 'identity') +
  theme_minimal() +
  labs(x = 'Assets', y = 'Weights', title = "Tangency Portfolio Weights") +
  scale_y_continuous(labels = scales::percent) 

ggplotly(p)

Then we can plot the efficient frontier and view how the portfolios fair on a complete scale

p <- portfolio_values %>%
  ggplot(aes(x = Risk, y = Return, color = SharpeRatio)) +
  geom_point() +
  theme_classic() +
  scale_y_continuous(labels = scales::percent) +
  scale_x_continuous(labels = scales::percent) +
  labs(x = 'Annualized Risk',
       y = 'Annualized Returns',
       title = "Portfolio Optimization & Efficient Frontier") +
  geom_point(aes(x = Risk,
                 y = Return), data = min_var, color = 'red') +
  geom_point(aes(x = Risk,
                 y = Return), data = max_sr, color = 'red')+
  annotate('text', x = max_sr$Risk, y = max_sr$Return, label = "Tangency Portfolio") +
  annotate('text', x = min_var$Risk + 0.025, y = min_var$Return, label = "Minimum variance portfolio") + 
  geom_point(aes(x = my_port_risk, y = my_port_returns), color = 'orange') +
  annotate('text', x = my_port_risk, y = my_port_returns, label = "My Portfolio")
 # once you find your port's risk and returns, can plot labels and point on graph 
 
ggplotly(p)

Thus we find the efficient frontier graph, be aware that this is only using past data and there is no prediction in this code.


Portfolio Analytics

Now that we have graphed the efficient frontier and seen how our portfolio compares to the optimized ones from 5,000 random trials. Let’s look at where we stand numerically, for some simplicity if we do not want to have to look at our random plots for each and every case.

We will look at the portfolio returns annualized, as well as it’s calculated risk, sharpe ratio, and how it compares to the benchmark. This may be more suitable for situations where we want to see how a single added holding influences these factors based on historical data.

Lets load our packages first

library(PortfolioAnalytics)
## Loading required package: foreach
library(fPortfolio)
## Loading required package: timeDate
## 
## Attaching package: 'timeDate'
## The following objects are masked from 'package:PerformanceAnalytics':
## 
##     kurtosis, skewness
## Loading required package: timeSeries
## 
## Attaching package: 'timeSeries'
## The following object is masked from 'package:zoo':
## 
##     time<-
## Loading required package: fBasics
## 
## Attaching package: 'fBasics'
## The following object is masked from 'package:TTR':
## 
##     volatility
## The following objects are masked from 'package:PerformanceAnalytics':
## 
##     kurtosis, skewness
## Loading required package: fAssets
# lets bind our portfolio 
portfolioPrices <- NULL
for (ticker in tick) {
  portfolioPrices <- cbind(portfolioPrices,
                           getSymbols.yahoo(ticker,
                                            from="2020-01-01", 
    periodicity = 'daily', auto.assign = FALSE)[,4])
  
  # us taking the close data
}

# check if missing data 

colSums(is.na(portfolioPrices)) # checking how many NA values over the time period
## SCHG.Close  RSP.Close SPLG.Close  IHI.Close  TSM.Close MSFT.Close  EFX.Close 
##          0          0          0          0          0          0          0 
##  FCX.Close BBJP.Close XLRE.Close  DIS.Close PYPL.Close AMZN.Close  PSA.Close 
##          0          0          0          0          0          0          0 
## GOOG.Close  VPU.Close 
##          0          0

Now we have our prices over the desired period, lets go ahead and find the returns. This gives us the portfolio returns each day as the ROC() function takes the rate of change at each close.

portfolioRets <- na.omit(ROC(portfolioPrices)) # want to omit the NA values

For our comparison to the benchmark we will be referring to the S&P 500, by using the SPY ticker.

benchmarkPrices <- getSymbols.yahoo('SPY', from = '2020-01-01',
                                    periodicity = 'daily',
                                    auto.assign = FALSE)
benchmarkRet <- na.omit(dailyReturn(benchmarkPrices$SPY.Close, type = "log"))

Now we have portfolio returns and the benchmark returns, now we want to aggregate the returns using our weights to show the actual port returns

portfolioReturn <- Return.portfolio(portfolioRets)
# this is using geometric returns, base is TRUE

Now that we have all this stuff, we can put together some of our metrics for the portfolio like beta and alpha. This is all pretty much using the performance analytics

print(CAPM.beta(portfolioReturn, benchmarkRet, .035/252))
## [1] 1.029876
# the risk free rate in this case is some number over 252
# 252 represents the number of trading days in a yr

print(CAPM.jensenAlpha(portfolioReturn, benchmarkRet, 0.035/252))
## [1] -0.02309981
SharpeRatio(portfolioReturn, 0.035/252)
##                               portfolio.returns
## StdDev Sharpe (Rf=0%, p=95%):       0.007061137
## VaR Sharpe (Rf=0%, p=95%):          0.004339316
## ES Sharpe (Rf=0%, p=95%):           0.002110014
# we want the STDDev Sharpe

print(SharpeRatio(portfolioReturn, 0.035/252)[1,])
## [1] 0.007061137

Now we can show some returns annualized and calendar year

table.AnnualizedReturns(portfolioReturn)
##                           portfolio.returns
## Annualized Return                    0.0309
## Annualized Std Dev                   0.2577
## Annualized Sharpe (Rf=0%)            0.1201
table.CalendarReturns(portfolioReturn)
##       Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
## 2020 -1.1 -0.4 -1.3 -1.0  0.8  1.7  0.3 -0.4  0.6 -1.5 -0.4  0.3
## 2021 -1.8 -0.7  0.8 -0.7  0.3 -0.3 -1.2  0.2 -1.0 -0.2 -1.9 -0.3
## 2022  2.3 -0.4 -1.6 -4.4 -0.3 -1.1  1.5 -0.5 -1.1 -0.9  3.8 -0.7
## 2023  1.5  0.2  1.6  0.5 -0.7  1.6   NA   NA   NA   NA   NA   NA
##      portfolio.returns
## 2020              -2.4
## 2021              -6.5
## 2022              -3.6
## 2023               4.9